function [val,x] = gurobi_milp(q)

    q = q'
    N = 6;
    M = 100000;

    % q = [1 1 1 1 1 1];
    % C = [0 12  M  M  9 16;
    %      12  0 19 12  M 15;
    %       M 19  0 21  M 17; 
    %       M 12 21  0 10 16;
    %       9  M  M 10  0 10;
    %      16 15 17 16 10  0];

    C = [0  7  7 10 10 12;  %Hardcoded N
         7  0  7 10 10 12;
         7  7  0 10 10 12; 
        10 10 10  0  7 12;
        10 10 10  7  0 12;
        12 12 12 12 12  0];

    c = [reshape(C,N^2,1);zeros(N^2,1)];% gurobi input: vectorized form of C Matrix

    objtype = 1;                        % gurobi input: 1 for minimize, -1 for maximize

    % Number of variables: integers N^2 for xij and N^2 binary for yij (could be N^2-N each)
    % Creating the R matrix required for one of the constraints

    R = (sum(q)-q(1))*ones(N,N);
    for i=1:N
        for j=1:N
            if(j==1)
                R(i,j) = q(1);
            end
            if(i==1)
                R(i,j) = sum(q);
            end
        end
    end


    A00 = [diag([1 zeros(1,N-1)]) zeros(N)];
    A01 = [zeros(N) diag([1 zeros(1,N-1)])];
    for i=1:N
     A0(i,:) = reshape(circshift(A00,[i-1,i-1]),2*N^2,1)'; 
     A0(i+N,:) = reshape(circshift(A01,[i-1,i-1]),2*N^2,1)';
    end

    A2(1,:) = [zeros(1,N^2) ones(1,N) zeros(1,N^2-N)];                  %colsum of yij
    A1(1,:) = [zeros(1,N^2) reshape([ones(N,1) zeros(N,N-1)]',1,N^2)];  %rowsum of yij
    for i=2:N
        A2(i,:) = circshift(A2(1,:)',N*(i-1))';
        A1(i,:) = circshift(A1(1,:)',i-1)';
    end

    A3 = [ones(N,1); zeros(N^2-N,1); zeros(N^2,1)]'; % the Nth leg flow value is 1 back to the first node. column sum


    A42(1,:) = [ones(1,N) zeros(1,N^2-N) zeros(1,N^2)];                 %colsum of xij
    A41(1,:) = [reshape([ones(N,1) zeros(N,N-1)]',1,N^2) zeros(1,N^2)]; %rowsum of xij
    A4(1,:) = A42(1,:) - A41(1,:);
    for i=2:N
        A42(i,:) = circshift(A42(1,:)',N*(i-1))';
        A41(i,:) = circshift(A41(1,:)',i-1)';
        A4(i,:) = A42(i,:)-A41(i,:);
    end
    bA4 = q';
    bA4(1) = bA4(1) - sum(q);

    A5 = zeros(N^2,2*N^2);
    for i=1:N
        for j=1:N
            A50 = [zeros(N) zeros(N)];
            A50(i,j) = 1;           % for xij
            A50(i,j+N) = -R(i,j);   % for yij
            A5(N*(i-1)+j,:) = reshape(A50,2*N^2,1)';
        end
    end

    randomA6 = [reshape([ones(1,N); zeros(N-1,N)],1,N^2) zeros(1,N^2)]; % the 1st leg flow value is N exiting from first node. row sum

    A =  sparse([A0; A1; A2; A3; A4; A5]);
    b = [zeros(2*N,1); ones(N,1);ones(N,1); q(1); bA4; zeros(N^2,1)];
    contypes = '===============================<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<';


    lb = zeros(2*N^2,1); % scalar means a uniform lower bound equal to scalar (which is zero here)
    ub = [sum(q)*ones(N^2,1);ones(N^2,1)]; % using loosely somewhat. Shoudl Rij figure here?
    vtypes = 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB';

    clear opts
    opts.IterationLimit = 4000;
    opts.FeasibilityTol = 1e-6;
    opts.IntFeasTol = 1e-5;
    opts.OptimalityTol = 1e-6;
    opts.LPMethod = 1;         % 0 - primal, 1 - dual
    opts.Presolve = -1;        % -1 - auto, 0 - no, 1 - conserv, 2 - aggressive
    opts.Display = 1;
    %opts.LogFile = 'weighted_fischetti_gurobi_mex_MIP.log';
    %opts.WriteToFile = 'weighted_fischetti_gurobi_mex_MIP.mps';

    [x,val,exitflag,output] = gurobi_mex(c,objtype,A,b,contypes,lb,ub,vtypes,opts);
    %if(exitflag==2)
    %    round(reshape(x(37:72),6,6))
    %end
    if(exitflag~=2)
        val = 10^6; % Arbitrary large number, might leade to errors if careless.
    end
end


% Syntax for Gurobi:
%     x = gurobi_mex(c, objtype, A, b, contypes, lb, ub, vartypes, options); 
%     *  c: objective coefficient vector, double. 
%     [] (empty array) means uniformly 0 coefficients, and scalar means all coefficients equal to scalar.  
%     * objtype: 1 (minimization) or -1 (maximization).
%     * A: constraint coefficient matrix, double, sparse.
%     * b: constraint right-hand side vector, double. 
%     * contypes: constraint types. Char array of '>', '<', '='. 
%     * lb: variable lower bound vector, double. 
%     * ub: variable upper bound vector, double. 
%     * vartypes: variable types. Char array of chars 'C', 'B', 'I', 'S', 'N'. C for continuous; B for binary; I for integer; S for semi-continuous; N for semi-integer. [] (empty array) means all variables are continuous. 
% Output Description
%     * x: primal solution vector; empty if Gurobi encounters errors or stops early (in this case, check output flag).
%     * val: optimal objective value; empty if Gurobi encounters errors or
%     stops early.